Mathematical analysis and simulations of the neural circuit for locomotion in lamprey 



Li Zhaoping^, Alex Lewis^, and Silvia Scarpetta^ 

^ University College, London, UK 
^ INFM & Dept. of Physics, University of Salerno, Italy 

We analyze the dynamics of the neural circuit of the lamprey central pattern generator (CPG) 
This analysis provides insights into how neural interactions form oscillators and enable spontaneous 
oscillations in a network of damped oscillators, which were not apparent in previous simulations or 
abstract phase oscillator models. We also show how the different behaviour regimes (characterized 
by phase and amplitude relationships between oscillators) of forward/backward swimming, and 
turning, can be controlled using the neural connection strengths and external inputs. 
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Locomotion in verbebrates (walking, swimming, etc.) 
is generated by central pattern generators (CPGs) in the 
spinal cord. The CPG for swimming in lamprey is one 
of the best known 0, Q , and has been a model system 
for investigations. It produces left-right anti-phase oscil- 
latory neural and motor activities propagating along a 
body composed of around 100 segments. A head-to-tail 
negative or positive oscillation phase gradient, of about 
1% of an oscillation cycle per segment, gives forward or 
background swimming respectively, and one wavelength 
from head to tail. External inputs from the brain stem 
switch the CPG between forward and backward swim- 
ming of various speeds and turning. Since isolated sec- 
tions of the spinal cord, down to 2-3 segments long 0, 
can produce swimming-like activity, the oscillations are 
thought to be generated by the neurons within the CPG. 
The neural circuit responsible is shown topologically in 
Fig. n It has ipsilaterally projecting excitatory (E) neu- 
rons and inhibitory (L) neurons, and contralaterally pro- 
jecting inhibitory (C) neurons, and provides output to 
motor neurons via the E neurons. All neurons project 
both intra- and inter-segmentally. The projection dis- 
tances are mainly within a few segments, especially from 
E and C neurons, and are longer, and possibly stronger, 
in the head-to-tail or descending direction 0-0 • 

Previous analytical work mainly treated the CPG 
as a chain of coupled phase oscillators in a general form 
()i — oji + J2j fij{(^ii ^j)- Here 6i is oscillation phase and 
LUi is intrinsic frequency, modelling the behaviour of one 
segment, and fij{9i, 9j) models inter-segmental coupling. 
This approach provided important insights into the con- 
ditions for phase-locked solutions applicable to various 
systems of coupled oscillators. However, its generality 
obscures the roles of specific neural types and their con- 
nections in generating and controlling behaviour. More 
recently, bifurcation analysis of the dynamics of a single 
segment was carried out, for a phase oscillator model de- 
rived from a kinetic (Hodgkin Huxley) equation for neu- 
rons 6] , and for a neural circuit model similar to the one 
used in this paper 0. Extensive simulations, including 
all neural types and detailed neural properties, have re- 
produced many features of experimental data |l| , though 
the model's complexity limits further understanding. 




FIG. 1: The lamprey CPG circuit. The solid and dashed lines 
denote intra- and inter-segment connections respectively. 

In all previous approaches, it is assumed that a sin- 
gle segment in the CPG can oscillate spontaneously, con- 
trary to experimental evidence that at least 2-3 segments 
are needed for oscillations We present an analytical 
study, confirmed by simulations, of a model of the CPG 
neural circuit in which an isolated single segment has a 
stable fixed point, with spontaneous oscillations occur- 
ring only in chains of coupled segments. The phase os- 
cillator approach is not applicable here since it assumes 
spontaneously oscillating individual segments perturbed 
by inter-segment coupling. Including specific cell types 
and their connections enables us to analyse the role of 
each of them in generating and controlling swimming. 
We show how external inputs select forward and back- 
ward swimming, by controlling the relative strengths of 
connections between various neurons, and produces turn- 
ing, by additional input to one side of the CPG only. We 
also analyse behaviour near the body ends. 

We model the CPG neural circuit with N — 100 
segments denoted by i = l,.,iV. The vector states 
(E;,Li,Ci) and (Er,Lj.,Cr), modelling the membrane 
potentials of the local populations of neurons at the 
left and right side of the body respectively, with E; = 
{E^i,E'^i,- ■ ■ ,E^i) etc., are modelled as leaky integra- 
tors of their inputs: 

E, = -E,-KV(a) + JV(EO+l£,i 

L, = -L,-A"5c(a)+WV(E0+lL,i (1) 

C, = -C, -BV(a) + QV(E/)-HV(Li)+Ic,i 

with the same equation for swapped subscripts r). 
.g£;(E;) = {gE{El), . . . , gE{E^)) are the neural activ- 
ities or firing rates, as non-negative (sigmoid-like) ac- 
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tivation functions of Ej, and likewise for gc{Ci) and 
5l(L,). KO, jo, AO, WO, B°, Q°, and are iV x 
matrices of non-negative elements modeling the synap- 
tic strengths between neurons. 1e,i, and Ic.i are 
external inputs, including those from the brain stem, as- 
sumed to be static. A left-right symmetric fixed point 
(E,L,C) where (E,L,C) = exists by setting external 
inputs to Ie,i = E; -I- K'^gc{Cr) - J"5£;(E;) (and analo- 
gously for other I's). Dynamics for small deviations from 
(E, L, C) can be approximated linearly, and, with a coor- 
dinate rotation (_E± , L± , C±)= [(E, , , C,) - (E, L, C)] ± 
[(Er, Lr, Cr) — (E, L, C)], transformed into two decoupled 
modes - the left-right synchronous mode (E_|_,L+,C+) 
and the antiphase mode (E_,L_,C_). Swimming re- 
quires oscillations, with wavelength of one body length, 
in the anti-phase mode while the synchronous mode is 
damped. The linearised equations are 



E± = -E± =F KC± 
L± = -L±tAC±- 
C± = -C±tBC± 



- JE± 
WE± 

QE± - HL± 



(2) 



where K =_K"g^(C), A = A°g'(j{C), B = B%(C), 
J = J05b(E), W = VJ^g'^m, Q = Q%(E), and 
H = H°g^(L) are effective connection matrices, and 
the g'(.)'s denote derivatives. The C neurons thus be- 
come effectively excitatory in the anti-phase mode. Not- 
ing that the lengths of the neural connections are much 
shorter than the body, and that isolated sections of 
spinal cord from any part of the body generate oscil- 
lations with similar amplitude and phase relationships 
P, , we make the approximation of translation invari- 
ance, so that matrix elements such as Jy depend only 
on (j — j), and impose the periodic boundary condi- 
tion, Jij = J(x), where x — {i — j) mod iV. This is 
adequate when behaviour near body ends is not con- 
sidered. Then all connection matrices commute with 
each other, with common eigenvectors (expressed as 
functions of segment number x) (E(a;), L(x), C(x)) oc 
g»(27rm/JV):r f^^. integer eigenmode -N/2 < m < N/2. 
The system solutions are thus combinations of modes 
(E±(a;,i),L±(x,t),C±(a;,i)) cx e^™*+»(2WAf)^ where 
is eigenvalue of eq. Q for mode m. Forward 
swimming results if the real part i?e(A^) < for all 
modes except the antiphase mode with m — 1, ie 
Re{Xi) > 0. Then this mode dominates the solution 
(whose growing amplitude will be constrained by non- 
linearity) (E(a;,t),L(a;,t),C(x,t)) cx e^-^^^Dt-^M-M, 
with oscillation frequency uj = |/m(Aj~)| and wave num- 
ber k — 2-k/N . Using the convention e^"^* for oscilla- 
tions, we omitted the solution cx e^'^'y^i )t+^^t jn the con- 
jugate pair of eigenvalues. To simplify our system, we 
note from experimental data that in forward swimming, 
E and L oscillate roughly in phase within a segment, 
while C leads them y]. We scale our variable defini- 
tions so that E_ = L_ in forward swimming. Then eq. 
imphes that (K - A)C^ = -(J - W)E_ in forward 
swimming. Since E and C have much shorter connections 



than the wavelength of oscillations during swimming, the 
connection matrices have zero elements far from the di- 
agonal, making (K — A)C_ and (J — W)E_ roughly either 
in phase or in anti-phase with C_ and E_ respectively. 
As C_ phase leads E_, (K - A)C_ = -(J - W)E_ is 
impossible unless (J - W)E_ = (K - A)C_ = 0. For 
simplicity we henceforth assume J = W and K = A, since 
non-swimming modes do not concern us. Consequently 
E± = L± and 



E±\ / J-1 tK 
C±j l-(H-Q) tB-1 



E± 
C± 



(3) 



where L and E are treated as a single population inhibit- 
ing or exciting C via connections H — Q. The eigenvalues 
for mode m are 



A™ 



-2 + Jm- Bm ± ^jRm + 2{Bl, + J2jj /2 

-2 + Jra + Bra - i/R^l /2. (4) 



Jm = Y.xKx)er''^'^''"'/^> is the eigenvalue of J (and 
analogously for other matrices), and i?,„ is the eigenvalue 
of R = 4K(H-Q)-(B-J)2. 

To elucidate the conditions needed for the antiphase 
mode with m = ±1 for forward or backward swimming 
to dominate, we analyse the bifurcations which occur 
as A^ for each mode (m, ±) changes as the effective 
neural connections are varied, either directly or via the 
external inputs. First, we focus on the left-right mode 
space (as in |^ Q for a single segment) of -I- and — , 
i.e., the synchronous and antiphase modes, by simply 
taking m — Q. Then, Jq, Bq, Hq, Kq, and Qq, are 
all real and non-negative, each being the total connec- 
tion strength on a postsynaptic cell from all cells of a 
particular type. Oscillation in the antiphase mode re- 
quires i?o > 0, necessitating > Qq, or that in the 
AC component of interactions above the background DC 
level, C neurons receive stronger inhibition from L neu- 
rons than excitation from E neurons. Consequently, 
Xq is real and the synchronous mode is non-oscillatory. 
As neural connections increase from zero, the antiphase 
mode undergoes a Hopf bifurcation when Re{XQ ) = 0, 
at Jo -l- = 2, and the synchronous mode undergoes a 
pitchfork bifurcation when Aj = 0, which occurs when 
(So -I- 1)(1 - Jo) = Ko{Ho - Qo). Oscillations result if 
the Hopf bifurcation has occurred but the pitchfork bi- 
furcation has not, i.e., Jie(Ag") > > X^ . The condition 
Aj < implies {Bq + 1)(1 - Jo) > Ko{Ho - Qo), ne- 
cessitating Jo < 1. Meanwhile, i?e(Ao ) > Ag leads to 
Bq > ^ Jq + Ro/2 > Jo, meaning that there must be 
sufficient inhibitory connections between left and right 
C cells. The J, W, and Q connections from E cells have 
to be relatively weak, consistent with the findings of . 
(If Rq < 0, the antiphase mode will undergo a pitchfork 
bifurcation, and the synchronous mode either a pitchfork 
or Hopf bifurcation. These regimes are less relevant to 
modeling the lamprey.) 
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Assuming the synchronous mode is damped, we focus 
now on the antiphase mode in the m mode space. Hopf 
bifurcations occur sequentially in various modes m in the 
order of descending i?e(A~). Taylor expanding J„i (and 
similarly _B,„, Rm) for small wave number k — 2Trm/N as 
is relevant for swimming, Jm = jo — ikji — fc^j2 + 0{k^) 
with jn =J2x-^(^)^' h^'^s 

2Re{X-{k)) = -2+jo+&o-fcri/(2V^)-fc2(j2+&2)+0(fc') 
2Im{X-{k)) = V^ + 0(fc) (5) 

making the mode with k ~ — ?'i/[4y^(j2 + &2)], which 
has the largest Re{X~{k)), dominant. From the defini- 
tion, (^0,^2, ^2) > 0, while stronger and/or longer con- 
nections in the descending direction imply > 0. 
Simply, J (and similarly for other matrices) is said to be 
descending (or ascending), if ji > (or ji < 0). Hence, 
if R is ascending, i.e., ri < 0, Re{X~{k)) increases with 
increasing fc, and the dominant wave number can be set 
to fc = 2tt/N for mode m = 1 by tuning the values of R, 
J, and B. If the connection strengths are such that only 
the TO = 1 mode undergoes the Hopf bifurcation, for- 
ward swimming emerges spontaneously. Switching R to 
descending leads to backward swimming. Note that J, B, 
H, K, and Q are all descending, multiplications and sum- 
mations of descending connections are still descending, 
and negating a descending connection makes it ascend- 
ing. Since B and H have to dominate J and Q respec- 
tively, R is composed of an ascending term — (B — J)^ and 
a descending term 4K(H — Q). Depending on the relative 
strengths of these two terms, R can be made ascending 
or descending to achieve forward or backward swimming. 
This could be achieved by changing the static inputs to 
shift the fixed point (E, L, C) of the system to a different 
gain regime (7^(E), (/^(L), (7^(C), and thus different ef- 
fective connection strengths H = H°(7^(L), etc. without 
changing the underlying connection structure H''. Al- 
ternatively, the external inputs might recruit extra func- 
tional cells to alter the effective connection strengths (sj . 

When connections are such that additional modes sat- 
isfy i?e(A^J > 0, the resulting behavior depends on the 
nonlinear coupling between modes. For illustration, con- 
sider nonlinearity only in gc{C)- 

E± = -e±tkV(c) + je± 
L± = -l±taV(c)+we± 

C± = -C±TB°.g±(C) + K'E±-HL± (6) 

where g±{C) ~ gc{Ci) ± gc{Cr)- If the nonlinearity is 
of the form gc{x) =2^-1- ax'^ — bx^ + 0{x'^), we have 

g_(C) kC-+ aC+C- - bCl/A - 3bC-Cl/4 (7) 
g+{C) « C+ + aCl/2 + - bCl/4 - ibC+Cl/A 

Hence, when C_ = 0, C+ cannot excite it since (?-(C) = 
0. However, if a ^ 0, the synchronous mode, will be 
excited passively by the antiphase mode through the 
quadratic coupling term aC^/2, responding with double 
frequency, as could be easily tested. 



To analyse coupling between the antiphase modes, we 
assume for simplicity that gc(C) is odd, so = 
since the synchronous mode is damped, and 5_(C_) = 
2gc{C- /2). Consider a small perturbation, in the to' 
mode direction, to the m = 1 cycle (the final orbit re- 
sulting from a small deviation from the fixed point in 
the TO = 1 mode, with a fundamental harmonic in the 
TO = 1 mode) such that C_(a;) k, Ci cos(27ra;/7V) -I- 
Cm' cos{2iTm'x/N) with Cm' ^ C*!- Expressing g-{C) 
as g^ (C) = f/„e'^''"'^/^, it can be shown that for large 
N , gm' ~ Cm'g'c where g'^ is the derivative of gc aver- 
aged over the unperturbed cycle. (More detailed analysis 
will be given in a future paper.) Because of the sigmoid 

form of gc{C), g'c < g'ciC). Then Cm' cx e^»'* with 
A^, as in equation Q except that {Bm',Km'), values 
derived from connections from C cells, are rescaled by 
a factor g'^/g'^{C) < 1. Thus the swimming cycle at 
large amplitude always remains stable against perturba- 
tion in other modes, even when the fixed point is unstable 
against these perturbations. 

If Re{Xi) ;» i?e(A~,) > 0, the to' cycle will have a 
small amplitude and hence g'{C) ~ g'{C) and it will be 
unstable against perturbation in the to = 1 mode. For 
larger i?e(A~/) the amplitude of the cycle is larger and ei- 
ther cycle will be stable. Suppose the neural connections 
are such that the to = ±1 cycles, giving forward or back- 
ward swimming, are both stable. The system would then 
display hysteresis, with the final behaviour depending on 
the initial conditions. Forward or backward swimming 
could then be selected by transient inputs from the brain 
stem, rather than by setting constant inputs as described 
above. This seems less likely to be the actual selection 
mechanism since experiments on Active swimming (pre- 
sumably with random initial conditions) seldom observe 
spontaneous backward swimming. However, the forward 
swimming could simply have a larger basin of attraction 
than backward swimming. 

When the lamprey turns, neural activities on left and 
right sides are unequal. This is realizable by adding an 
additional constant input to one side in the animal and in 
models (Sj, |S] , leading to unequal mean activities without 
disrupting the oscillations, provided that the gains g'{.) 
are roughly constant near the fixed points. Simulations 
(Fig. 2) confirm the analysis above. 

To study behaviour at the body ends or in short sec- 
tions of spinal cord fl|, or equivalently to see the effects 
of longer connections, we abandon translation invariance. 
Eliminating C in eq. ||2Jl, the minus mode has: 

E + (2 - J - B)E + [1 - J - B + BJ + K(H - Q)]E = 0, 

or, oscillator i is driven by force Fi from other oscillators 

E, + (2 - J,, - B,,)E, + [1 - R,,]E, 

where R = B+J — BJ — K(H — Q). The intrinsic oscillation. 
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FIG. 2: Simulations. A: Membrane potentials of E popula- 
tion on either side of one segment during forward swimming 
and turning. The oscillations are in anti-phase between the 
two sides. Turning is induced by an additional constant input 
to one side only, starting at the time indicated by the dashed 
line. B: C slightly phase leads E during forward swimming. C 
& D: Waveform of E along the body in forward and backward 
swimming, at consecutive times increasing in the direction in- 
dicated by the arrows, in the translational invariant model. 
The switch from forward to backward swimming is achieved 
by increasing the strength of H and Q. E: Oscillation wave- 
forms (note different amplitudes) in body segments at head, 
tail and centre of the body, without translational invariance. 
F: Forward and backward swimming from different initial con- 
ditions, with the same connection strengths and inputs. 



Ei ~ is damped, i?e(A) = -1 (J + 6)^^/2 < 0, as 
indicated by experiments 0, Q . We estimate Fi using 
the approximation that oscillators j ^ i still behave as 
Ej cx e-'('^*-'=-''). We then have Fi = a.'Ei + /J^E,, where 



= (J + cos(/c(j - i)) - Rij sin(fc(j - i))/uj 

Pi = J2j^i (J + B)ytjsin(fc(j - i)) + Rij cos{k{j - i)) 



The term a^Ei when ai > feeds oscillation energy 
into the i*'' (receiving) oscillator, causing emergent os- 
cillations in coupled damped oscillators. We divide 
~ ctj.dosc + cti.asc into the descending and ascend- 
ing parts, with summations over J2j<i ^^'^ ^j>i 
spectively. Hence, for i = I, ai — a^^asc; for i = N, 
ctN = ai,dcsc7 and for 1 <^ i <^ N , ai = ai + ■ Since 
the first and last segments oscillate due to the driving 
force from other oscillators, ai > and > 0. Conse- 
quently, ai < a-^ji and ajsi < cx]sj/2- Further, since de- 
scending connections are stronger, it is most likely that, 
for 1 ^ i < iV, ai,desc > oti,asc- Consequently, ax < un- 
Hence, the rostral oscillator has a smaller amplitude than 
the caudal one, which in turn has a smaller amplitude 
than the central one (Fig 2(E)). Firing rate saturation 
and variations of the fixed point along the body may ob- 
scure this pattern in experimental data, although body 
movements are indeed smallest near the head |T]|. Sim- 
ilarly, oscillation amplitudes will be reduced in sections 
of spinal cords shorter than the typical lengths of inter- 
segment connections, and will eventually be zero in ever 
shorter sections, as observed in experiments 0- 



In summary, analysis of a model of the CPG neural 
circuit in lampreys has given new insights into the neu- 
ral connection structures needed to generate and con- 
trol the swimming behaviour. In particular, we predict 
that the contra-lateral connections between C must be 
stronger than the self-excitatory connection strength of 
the E neurons; that the C neurons are more inhibited 
(in their AC components) by L neurons than excited by 
the E neurons; and have shown how different swimming 
regimes can be selected by scaling the strengths of the 
various neural connections without changing the connec- 
tion patterns. Our framework should help to provide 
further insights into CPGs of animal locomotion. 
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